%% =========================== [Description] ============================
%% ---------------- implements the function -----------------
% Under single input
% Draw the mean variance curve of acceleration displacement         
% Draw a probability density curve for a given second             
% Draw the response comparison curve for a deterministic and random input
% ---------------- predecessor file -----------------
% quasi0StiffSuspensionCertain.slx
% quasi0StiffSuspension.slx
% n1PolyRoots.m
% GetStatistics.m
% ------------------ enter -------------------
% Cylinder air pressure
% Piston distance from rear wall
%% ==================== [Select excitation signal and output] =====================
clear
disp('Pavement sine signal - input 1')
disp('Pavement random signal - input 2')
exc=input('Signal:');
disp('Output select acceleration - input 1')
disp('Output Select Displacement - input 2')
outputAmount=input('outputAmount:');
swiMat=[-1,1];
excSwitch=swiMat(exc);
%% ===================== [model run to find distribution] ==========================
%% ----------------------------- [Determine values] ----------------------------------
d0Certain=32.5e-3; % piston distance from back wall m
mCertain=800; % mass on the spring kg
AcCertain=pi*0.1^2; % area of cylinder piston m^2
L0Certain=0.467; % cylinder node distance m
kCertain=3; %cylinder air pressure atm
m1=100; %Tire mass
%--------------------------------------------------------------------------
variableCell={d0Certain,AcCertain,L0Certain,kCertain};
disp('d0=1,Ac=2,L0=3,k=4')
varIndex=input('Random variable:');
inVar=variableCell{varIndex};                                
delta=inVar*0.2; %standard deviation of input variables
order=3; % number of chaotic polynomials
Collocation=n1PolyRoots(order); %Find the matching points
iinVar=Collocation*delta+inVar; %standard normal distribution to relate to input variables
variableCell{varIndex}=iinVar;
m=800*ones(length(iinVar),1);

%------------------------------ [parameter] ------------------------------------
d0=variableCell{1};
Ac=variableCell{2};
L0=variableCell{3};
k=variableCell{4};
%--------------------------------------------------------------------------
Pa=101*10^3;
A_level=16; %-|
B_level=64; %-|-pavement
C_level=256; %-|
Level=B_level;
v=60; %Vehicle speed
C=2000; %damping
%--------------------------------------------------------------------------
sim('quasi0StiffSuspensionCertain.slx')
sim('quasi0StiffSuspension.slx')

switch outputAmount%Select the output amount
    case 1 
        outputData=acc.data;
        opDataCertain=accCertain;
        titleText='acceleration';
        yLabelText='Acceleration (m/s^2)';
    case 2
        outputData=dev.data;
        opDataCertain=devCertain;
        titleText='Displacement';
        yLabelText='displacement(m)';
    otherwise
        error('Select 1 or 2');
end
pointNumber=100;
timeNumber=10;
step=timeNumber/pointNumber;
time=0:step:timeNumber-step;
MCnumber=100;
PCmean=zeros(pointNumber,1);
PCvar=zeros(pointNumber,1);
MCdata=zeros(pointNumber,MCnumber);
MCinput=zeros(pointNumber,MCnumber);
for i=1:pointNumber    
    [PCmean(i),PCvar(i),~,~,MCdata(i,:),MCinput(i,:)]=GetStatistics ...
    (Collocation,outputData(100*i-99,:)',order,MCnumber); 
    %----------------------------------------------------------------------
    proPers=floor(i/pointNumber*100);                                 %
    clc                                                               %---progress bar
    disp(['[','#'*ones(1,proPers),' '*ones(1,100-proPers) ...         %
        ,']',num2str(proPers),'%'])                                   %
    % ----------------------------------------------------------------------
end
%% ====================== [Standard deviation and mean plot] =========================
figure(1)
errorbar(time,PCmean,sqrt(PCvar),'-o',...
    'color',[0,0,0],...
    'linewidth',0.8,...                        % Set line width     
    'MarkerSize',4,...                         % Set marker symbol size
    'MarkerEdgeColor','red',...                % Set marker edge size
    'MarkerFaceColor','red') %set the marker inner size
title([titleText,'Mean and standard deviation over time'])
xlabel('Time(s)')
ylabel(yLabelText)
% %% ======================[timeline probability density] ============================
% Xpdf=zeros(pointNumber,100);
% X=zeros(pointNumber,100);
% for k=1:pointNumber
% [Xpdf(k,:),X(k,:)] = ksdensity(MCdata(k,:)); 
% end
% for u=1:100
% for t=1:pointNumber
% if Xpdf(t,u)>30
% Xpdf(t,u)=30;
% end
% end
% end
% figure(2)
% pcolor((0:5)'*ones(1,5),ones(6,1)*(-0.2:0.1:0.2),zeros(6,5))
% hold on
% pcolor(time'*ones(1,100),X,Xpdf);
% shading interp; 
% colorbar; colormap(jet);
% xlabel('X');ylabel('Y');
%% ======================== [probability distribution for a given second] =========================
t=5;
tNumber=t/step+1;
[Xtpdf,Xt] = ksdensity(MCdata(tNumber,:));
figure(3)
plot(Xt,Xtpdf,'LineWidth',1)
title(['output',num2str(t),'probability density at s'])
xlabel('Output random variable')
ylabel('probability density \it{f}')
%% ==================== [chaotic mean vs model mean plot line] ====================
figure(4)
plot(time,PCmean,'-o','MarkerSize',4)
hold on
plot(opDataCertain.time,opDataCertain.data,'-.' ,'MarkerSize',4)
title([titleText,'Curve over time'])
ylabel('Output random variable mean')
xlabel('time(s)')
legend('random input','deterministic input')





















